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^ ■ General orbits of a particle of small mass fi around a Kerr black hole of mass M are 

' characterized by three parameters: the energy, the angular momentum and the Carter con- 

1/^ I stant. The time-averaged rates of change of the energy and the angular momentum can 

, be obtained by computing the corresponding fluxes of gravitational waves emitted by the 

fSJ ' particle. By contrast, the time-averaged rate of change of the Carter constant cannot be 

' expressed as a flux of gravitational waves. Recently a method to compute this rate of change 

, was proposed by Mino, and we refined it into a simplified form. In this paper we further 

. extend our previous work to give a new formulation without the aid of expansion in terms 

"'^ ' of a small inclination angle. 

7'. 

5--I . SI. Introduction 

Thanks to recent advances in technology, an era of gravitational wave astron- 
^ . omy has almost arrived. There are already several large-scale laser interferometric 

H I gravitational wave detectors in operation. Among them are TAMA300,i) LIG0,2) 

GEO-6003) and VIRGO. The primary targets of these ground-based detectors are 
inspiralling compact binaries, which are expected to be detected in the near future. 

There are also projects for space-based interferometric detectors. LISA is now 
on its research and development stage, and there is a future plan called DE- 
CIGO/BB0.6)'7) These space-based detectors can detect gravitational waves from 
solar-mass compact objects orbiting supermassive black holes. To extract physical 
information concerning such binary systems, it is essential to know the theoretical 
gravitational waveforms with sufficient accuracy. The black hole perturbation ap- 
proach is most suited for this purpose. In this approach, one considers gravitational 
waves emitted by a point particle that represents a compact object orbiting a black 
hole, assuming the mass of the particle is much less than that of the black hole 
M, i.e. // < M. 

To lowest order in the mass ratio /u/M, the orbit of the particle is a geodesic on 
the background geometry of a black hole. Already at this lowest order, combined 
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with the assumption of energy and angular momentum balance between the emitted 
gravitational waves and the orbital motion, this approach has proved to be very 
powerful in evaluating general rclativistic corrections to the gravitational waveforms, 
even for neutron star-neutron star binaries, for which the assumption of this approach 
is maximally violated.^) 

However, the deviation from the geodesic cannot be completely specified by the 
rates of change of the energy and angular momentum. In order to describe general 
orbits, we need to know the evolution of the third "constant" of motion, i.e. the 
Carter constant Q. For this purpose, we need to evaluate the gravitational self- force 
acting on the particle directly. Here, the gravitational self-force is the force due to 
the metric perturbation caused by the particle itself. 

Roughly speaking, there are two levels in the computation of the self-force in 
the linear approximation.^-* The advanced level is the direct computation of the 
time-dependent self-force, in which one computes the self-force without any further 
approximation. The main problem here is that the full (bare) metric perturbation 
diverges at the location of the particle, which is assumed to be point-like, and hence 
so does the self-force. About a decade ago, a formal expression for the gravitational 
self-force was found, which contains the expression for the Green function divided 
into two parts; the direct part and the tail part."'^''^' Later, this formula was 
reformulated in a more sophisticated manner by Detweiler and Whiting. In this 
new formulation, the direct part is replaced with the S part and the tail part with 
the R part. The R part has the improved property that it is a solution of source-free 
linearized Einstein equations. Thus, what we have to do to obtain a meaningful 
self-force is to compute the R part of the metric perturbation. The "mode sum" 
scheme, a practical calculation method for the R part, has been developed. 
There are several implementations of this scheme in a scalar toy model: Ref. 20) for 
a static particle, Refs. 22) and 21) for radial infall, and Refs. 23), 24), 25) and 26) 
for a circular orbit. Extensions to the gravitational case are made in Ref. 27) for 
radial infall, and in Ref. 28) for a circular orbit. 

Despite recent progress in the direct computation of the self-force, there seem 
to remain obstacles in computing it for general orbits on the Kerr background. Fur- 
thermore, all results obtained from the computation of the self-force are not equally 
meaningful. A few time-averaged combinations composed of the self-force are impor- 
tant for the prediction of the gravitational waveform. The most important quantities 
are the time-averaged rates of change of the "constants" of motion. Here a constant 
of motion means the quantity, such as the energy or angular momentum, that stays 
constant for background geodesies. Hence, the other, easier level of computing the 
self-force is to compute these rates of change of the "constants" of motion. Many 
years ago, Gal'tsov^^^ advocated using the radiative part of the metric perturbation 
to calculate dE/dt and dL/dt. The radiative field is defined as half the retarded 
field minus half the advanced field. The advantage of using the radiative field is that 
this force is relatively easy to compute and is free from divergence. Hence we do not 
have to worry about how to subtract out the divergent part. Gal'tsov showed that 
the calculation using the radiative field correctly reproduces the results obtained by 
using the balance argument for dE/dt and dL/dt when they are averaged over an 
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infinitely long time interval. Recently, Mino directly verified the validity of applying 
the same scheme to the computation of the time averaged rates of change of the 
constants of motion including the Carter constant.^'') 

From the equivalence of the two calculations for dE/dt and dL/dt shown by 
Gal'tsov, we see that, in order to obtained time-averaged rates of change of the 
energy and angular momentum, all we have to do is to compute the fluxes of the 
energy and angular momentum evaluated at infinity and on the black hole horizon 
by using the Teukolsky formalism. Therefore, it is expected that the same might 
be true for the Carter constant. Recently, we have shown that this is indeed the 
case, obtaining a new simplified formula of the time-averaged rate of change for 
dQ/dt written in terms of the asymptotic amplitude of gravitational waves. '^^^ This 
formula requires some knowledge of the particle orbit. In this sense, the rate of 
change is not expressed as a flux determined by the asymptotic form of gravitational 
waves. In our previous paper, using this new formula, we gave analytic expressions 
for the time-averaged rates of change of the energy, the angular momentum and the 
Carter constant in the post-Newtonian expansion. '^^^ There, we also made use of the 
expansion in terms of the orbital inclination angle for a technical reason. Here, we 
extend our previous results, eliminating the limitation to small inclination angles. 

This paper is organized as follows. In ^we review the basic formalism about 
how to analytically compute the rates of change of the constants of motion due 
to gravitational wave emission for a particle orbiting a Kerr black hole, which is 
developed in Ref. 32). The formulas for the rates of change are written in terms 
of the amplitude of each partial wave of the emitted gravitational waves. In ^we 
explain our new method for the analytic evaluation of this amplitude for a general 
orbit with a large inclination angle. In ^ substituting the results obtained in the 
preceding section into the formulas described in ^ we compute the time-averaged 
rates of change, dE/dt, dL/dt, and dQ/dt. We also compute the phase evolution of 
gravitational waves. In ^we summarize this paper. 

§2. Basic formulation for adiabatic radiation reaction 

In this section we give a brief review on the Teukolsky formalism'^^^'^'') as well as 
the basic formulas obtained in Ref. 31). We consider the background Kerr spacetime 
in the Boyer-Lindquist coordinates: 



ds^ = - il — j dt^ — dtdif + —dr^ 

+Ude^ + (^r^ + a^ + sin2 0^ sin^ 9d<p\ (2-1) 

where 

X: = r^ + a^cos^9, A = r'^ - 2Mr + . (2-2) 



Here, M and a are the mass and the angular momentum of the black hole, respec- 
tively. In the Teukolsky formalism, the gravitational perturbation of a Kerr black 
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hole is described by a master variable ^, which satisfies the master equation 



(2-3) 



where 



A 



2 ■ 2 

a sm I 



1 



^2 4Mar 



A 



A 



sm 



+A-'dr{A'+'dr) + —-de{sinede) + 2s 
sm d 

M(r^ - a2) 

- — r — la cos 



a{r — M) i cos ( 



A 



A sm" 

2 , 



dt - s{s coVe- 1), 



(2-4) 



and sT is the source term. The master variable ^ is equal to ipQ for s = 2 and 
(r — mcos0)^V'4 foi' = —2, where 



(2-5) 



are the so-called Newman-Penrose quantities. Here, Cap-yS is the Weyl tensor, and 
the null vectors i,n,m are defined by 



rri' 



t = ((r^ + a^),A, 0, a) /A, = ((r^ + a^), -Z\, 0, a) /{2U), 

= (msin6l,0,l,i/sin6l) / {V2{r + ia cos 6)) . (2-6) 



The bar denotes complex conjugation. 

The master equation (12-3P can be solved by decomposing the master variable ip 

as 



A J 



-iut 



(2-7) 



where A represents a set of separation constants, rn-,^}. The equations for the 
radial and angular parts can be separated, and we obtain 



A 



d_ 
dr 
1 d 

^[ede 



d 



+ 



A'+^ — 
dr 

9± 
dO^ 

m + s cos^ 6*)^ 



- 2is{r - M)K 



A 



+ Aisujr — A 



sin( 



2 2-2 

a uj sm I 



2aujs cos 6 + s + 2amoj + A 



sm 



,i?yl(r) = sTa, 



,Sa{0) = 0, (2-^ 



where 



K = (r^ + a^)a; — ma, 



and A is the eigenvalue determined by the equation for sSa- The angular function 
sSa is called the spin- weighted spheroidal harmonic, which is usually normalized as 



sSa) sinO d9 = 1. 



(2-9) 
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We define two independent homogeneous solutions of tlie radial Teukolsky equa- 
tion: 



, B'r''A-'e-'^'\ for r^r+, 

nup \ e + ZA e , ror r ^ r+, 

s^A I ^-2.-l^tranSgic.r*^ for r ^ +00. ^ ^ 

Here, k = u — ma/2Mr+, and r* is the tortoise coordinate defined by 

='^ + —In In— r— , (2-11) 



with r± = M ± V — a"^. A systematic analytic method to compute homogeneous 
solutions for sRa and sSa in the post-Newtonian expansion has been developed in 
Refs. 35) and 36). 

The explicit form of the source term Ty[ is given by 

Ta = 4 f dQdtp-^p-^{B'^ + B'*)e-''^'^+^^=^, (2-12) 
J y2iT 



where 



-^^p^-pA^L., [p-'fMp-'-p-'^-'T^n)] , 
B'2* = --/liA^J^ [p-^J+{p-''pTrrmi)\ 

^ /pZ\2j+ [p-V/A-i£_i(p-V'7mn)] , (2-13) 



2V27r 
with 

p = (r — ia cos)~^ , Cg = dg + — — - — awsin^ + scot^, J+ = dr -\ — t-. 

smO A 

In the above, T„„, 7m„ and are the tetrad components of the energy-momentum 
tensor. Here and hereafter, for simplicity, we restrict our attention to the s = —2 
case, and we omit the subscript s. 

We solve the radial Teukolsky equation using the Green function method. A 
solution of the Teukolsky equation that is purely out-going at infinity and is purely 
in-going on the horizon is given by 

Ra = ^S^Rj£ dr'B'j^TAA-' + R^I^"dr'RjTAA--'Y (2-14) 

where the Wronskian Wa is given by 

Wa = 2iuj&X^'^Bt''. (2-15) 
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Then, the solution has the following asymptotic behavior near the horizon: 

ptrans /1 2 „— jfer* poo 

""'^^ ^ ^ luCpBT L " Z^A^-^^r^. (2.16) 

Here we defined , the amplitude of each partial wave labelled by A. In the r — oo 
limit, the solution takes the form 

R^^r ^ oo) ^ / dr':^^ = Zfr^e"-'*. (2-17) 



Hereafter, we focus on the gravitational waves emitted to infinity, partly because 
the horizon in-going wave does not give the dominant contribution in the present 
calculation, and partly because the extension is rather straightforward. 

We consider the motion of a point particle, whose coordinates are = {tzir), 
i^z{t), Oz{t), ifzir)). Here r is the proper time along the orbit. For geodesic motion, 
there are three constants of motion. 



2Mrz \ t 2Marz sin^ 6z 



-u 



C = u""^^' = — + sm^ 9zu'^, 

sin 9z 

where u"' = dz'^/dr is the four velocity, and K^/s is the Killing tensor defined by 

Kap = 2Ue^^np)+r^ga0. (2-19) 

We often use alternative notation for the Carter constant, C = Q — {a£ — C)'^. For 
orbits on the equatorial plane, C vanishes. Using these constants of motion and a 
new parameter along the trajectory A defined by dX = dr/E, the geodesic equations 
become 

'^"'^ ^ ' R{rz), ^ = -a{aS sin^ 9z-C) + !i±^P(rz), 



where 



dX J ' dX ' ^ ' ■ A 

2 

d cos 6*2 \ dipz La 



0{cos9z), ^ = -aS + -^ + -P{rz), (2-20) 



P(r) = £{r'^ + a^) - aC, R{r) = [P{r)f - A\r^ + Q], 

©(cos 9)=C-{C^a^{\- E^) + C?) cos^ 9^a^(\- E^) cos^ 9. (2-21) 

It should be noted that the equations for the r-component and the ^-component are 
decoupled when they are written in terms of A. Both R and Q are quartic functions 
of their arguments. Hence, both solutions are given by elliptic functions. Taking 
the amplitude of the radial oscillation as a small parameter, we can systematically 



Adiabatic Evolution of Three 'Constants' of Motion for Greatly Inclined Orbits 7 



expand the radial solution r^^X) as a Fourier series. For the motion in the ^-direction, 
perturbative solutions in Fourier series can be systematically obtained even for a large 
inclination angle by taking the coefficient of the quartic term a^(l — i?^) as a small 
parameter. As we explain below, this expansion is a part of the post-Newtonian 
expansion. We stress that it is not necessary to restrict the amplitude of oscillation, 
\9 — tt/2\, to small values in this expansion. 

The other two equations in (l2-20p are integrated as 



' dip, 
dX 



(^,(A)=9JW(A)+(^W(A) + 



(9)1 



(2-22) 



where (• • • )a = lim^A^oo(2/iA) ^ / 



1 r-^-^ 



AX 



dX- 



represents the time average along the 



geodesic, and t^^^ and t^^\ which are defined by 

(r2 + a2)P(r,) / {r^, + a^)P{r,) 



t^''\X) = I dX 



A 



A 



t(^)(A) = - J dX{a^£ sin'^e, - {a^£ sin^ 9,) ^} 



(2-23) 



are periodic functions of periods 27r/i7r and 2Tr/f2g, respectively. The functions 99*-''-' 
and (f^^^ are defined in the same way. 



^M(A) 



V^(^)(A) = / dX 



dX 



aP{r,) 


/aP{r,) 


A 


\ ^ 


C 


1 \ 


sin^ 9z 





(2-24) 



Using Eqs. (j2-20p . the tetrad components of the energy momentum tensor are 
expressed as 



fi^6{r-r{t))5{9-9{t))5{ip-ipit)), 



'^mn l^ 



smt 

Crni 

sin ( 



5(^r-r{t))5{9-9{t))5{^-^{t)), 



^^S{r-r{t))6{9-9{t))6{^-^{t)\ 
sm 9 



(2-25) 



where 

Cmn 
Cmm 



At 



dr 



n 2 



2^/2i 



8{r^ + a") - 



dr 
dX 



isni9 ( a8 



C 



sm 



1 d cos ( 



sm ( 



dX 



P_ 
2i 



ism 9 I a£ 



sm 



1 d cos ( 
sin 9 dX 



(2-26) 
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and t = dt/dX. After some calculation, Eq. (j2-12p finally becomes 

TA = f^J dte'^'-'^'P^'^A' [[Ann^ + ArnnO + Arnrn^)5[r - r{t)) 



where 
and 



+ {{Amni + Afnfni)5{r - r(t))}_^ + {(^^^2)<5(r - r{t))} . 
C}^ = de : — - + auj sin 9 + a cot 6, 



(2-27) 



sin( 



A-nnO 
AfanO 

A-rnfhO 
A-ffinl 
■Aififhl 
■Afhfh2 



c- 



1 



{^\Sa) [—T+P + P]- as\neSA-r{p - p) 



P ^ pCmmSA 



A 
K 



K 



A 



A 



K 



Z\2 



A 



V^A 



p-^C^ 



P ^pCmfhSA P 



C^Sa + iasm9{p — p)Sa 



K 



A 



P pCmmSA- 



(2-28) 



Inserting Eq. into Eq. (pTT]) . we obtain Za as 



dR']^ 
dr 



mnl + ^mml } + 



d^R^ 
dr"^ 



-A 



mm2 



(2-29) 



r=r{t) 



Furthermore, when we consider bound orbits, the frequency spectrum of Ta becomes 
discrete. Therefore Za takes the form 



Za= 2vr<5(a; - u;^'''^)Zj^Vi^uj, 



(2-30) 



nr,ng 



with 



nr,ng 
m 



dX 



m 



d<Pz 
dX 



V 



(2-31) 



Here A denotes the set of parameters {i,m,nr,ng}. Then the time-averaged lumi- 
nosity is given by 



d£_ 
lit 



A 



(2-32) 
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and the time-averaged angular momentum flux is given by 



m 



(2-33) 



Here (• • • ) represents the time average with respect to the coordinate time t. We 
quote the formulas for calculating (dC/dt) from our previous paper'^^) : 



d£ 
lit 
dQ 
dt 



dQ 
dt 

r^ + a^)P 



2{a£ -£)\^a 
d£ 



d£_ 
It 



dC 



A 



dt 



A /^\dt 



nr,ng 



(2-34) 



§3. Amplitude of partial wave 



As mentioned above, taking the amplitude of the radial oscillation and the coeffi- 
cient of the quartic term of ©(cos 9), which are denoted by e and eo = a^(l — 
respectively, as small parameters, we can systematically expand the solutions as 
Fourier series. Indeed, the solutions are found to take the forms 



rJX) =p[l + e cos nrX]+0(e'^), cos 6^ JA) = J sin I?eA + ©(eo), (3-1) 

V y + 1 

where p and y are defined in Eqs. p-2p and (j3-3p . There is a difficulty when we ac- 
tually calculate the amplitude of a partial wave using Eq. (j2-29p . Here we want 
to perform the inverse Fourier transformation analytically. This is possible when 
Zy[ is expanded as a power series in sinusoidal functions. However, and 
have factors of 1/ sin^^ and 1/ sin^ 6^ through Cmn and Cmm, respectively. To make 
matters worse, e*"^ is proportional to l/sin0. This means that Zj^ has a term pro- 
portional to Y^(?/ + !)/(! + y cos^ apparently. It would seem that this would 
prevent us from performing the inverse Fourier transformation analytically, unless 
we expanded this factor with respect to y, assuming a small inclination angle. If such 
problematic terms remain, however, it follows that we cannot truncate the series ex- 
pansion labelled by ng at a finite order for general orbits. This seems quite counter 
intuitive. We expect, therefore, that we can somehow overcome this difficulty. In 
this subsection, we rewrite the source term into a tractable form, which does not 
contain any inverse power of sin^. We present the results up through 0(e^,f^) as a 
demonstration, where u is a typical velocity of the particle defined more precisely in 
Eq. 

3.1. General geodesic orbits in Kerr spacetime 

Generic Kerr geodesies can be specified by fixing three orbital elements, say, the 
semi-latus rectum p, the eccentricity e, and a dimensionless inclination parameter y. 
We define p and e such that the turning points of the radial motion, the apastron 
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and periastron, are 



P P lo o\ 



respectively. We also define y by 

y = c/c'', (3-3) 

which is something like the squared tangent of the inclination angle. This parametriza- 
tion of orbits is useful for obtaining an intuitive understanding, but it is, of course, 
equivalent to specifying the three "constants" of motion £, C and C. Solving the 
above defining equations, £ and C can be expressed in terms of p, e and y as 

f = 1 _ 1 2 + 3 4 _ ^^^5 ^ ^2 r 1 2 _ 3 4 ^ 

2 8 [^2 4 

P 2 \ 8 J 

■^ ,3 ^\^2„,4 , / 2-v^3 , 9^ V,5 



where 
and 



V = v/M/^, (3-5) 

y = / = , . (3-6) 

Then, expanding the periodic parts with the period of the radial oscillations in the 
geodesic equations in powers of v and e, the solutions up through 0(e^, f ^) are found 
to be given by 

'T-max 

r(A) =p COSrirf^r^j 

nr=0 

^max 

tW(A)= J^tWsinn.r^A 

nr=0 

^max 

^M(A)= ^ ^Wsinn,r2,A, (3-7) 

nr=0 

where nmax = 2, which is the truncated order in e, and 



ai = e 



1 i,2^„x.„,3 (i-2n9^A.,4 



ao = l + e^ -v^ + qYv" - ( 3 - ^ ^ ) + lOgy^^ 



Iq — u, 



2 2 
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= e {2 + - 6 qYv^ + (l7 - (1 - AY^) q^) - 54 qVv^ } , 



4 4 4 V 8 8 / 4 

(^J'') = e { -2^?;^ + 2 q^Yv"^ - 10 } , 

4')=.^|-2^„« + p}. (3.8) 

The angular velocity fir is obtained simultaneously when we solve r(A), but we give 
it below, together with the angular velocity of the ^-oscillations. 

The other parts of the geodesic equations which are periodic with the period 
of the 0-oscillations are expanded in powers of eg = a? {1 — E"^) / C"^ = 0{v^) and 
obtained as 

^max 

cos 0(A) = a/i — y2 ^ fj^^ smnof2gX, 

"9=0 

t^'HX) = E sinn,r?,A, (3-9) 

ng=0 

where rimax = 2x1 + 2 = 4, which is twice the truncated order in eo plus two, 

/3o = 0, Pi = l + eo^- 132 = 0, (3s = eo^——^, ^4 = 0, 

lb lb 

ff^o. .f'^^o^^^ii^. (3.10, 

For (f^^^ , we introduce the variable X = sin 6 e*'^**' . Then, it can be expanded as 

3?(X) = ^ Xl cosngneX, = ^ X^ smnef2gX, (3-11) 

where 

v« (^ + 1) , 1-^(1 -9y)(i-n ^ 

^ (l-Y) YHl-Y') ^ ^ YHl + Y){l-Yr 



Using Eq. (|3-4p . we re-expand the above solutions in powers of t; and e to obtain 



^max "^max 



cos 0(A) = Vl-y2 ^ sinn^I^eA, t(^)(A) = ^ t^^] smnefieK 

ng=0 ng=0 
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TT-max 



3fi(X) = J2 ^ne cos nef2gX, Q{X) = ^ smn0f2gX, (3-13) 

ng=0 ng=0 



where 

Po = 0, Pi = l + 

P3 = 



Ad) 



-^0 



(1 - 9Y^)<f 4 2 r (1 - 9Y^)<f 4 



-V — e 



16 t 16 

(1 - y2)^2 ^ ^ (1 _ Y2jq2 



P2 = 0, 



16 



— e 



16 



/?4 = 0, 



0, ^1 ^ — 0, ■ — ^ ~ 

+ (i-9y)(i-yV 

2 32 



(1 - Y')q' 3 (1 - Y')q^ , <e) 



v\ t'"' = 0, t 



(9) 



X. 



(i-9y)(i-y^)g^ 

32 

y(i-y2)g2 



Xf = 0, 



3? 



(i+y)(i-F)v 4 2r(i+n(i-nv 4 



32 



-V — e 



32 



c,_ (1-y) ^ (i + 5y)(i-yV ^^.4 ^2 



Xg^ = 0, 



16 



(i + 5y)(i-yv 4 



16 



(i + y)(i-y)V 4 , 2r(i + >^)(i-nV 4 



32 



-V + e 



32 



(3-14) 



Finally, the frequencies of radial and ^-oscillations and the non-oscillating parts 
of dt/dX and dip/dX are given by 



f2r 



1 / dt 



1 /dip 



p2 \dX/. M 



1- -v"^ + 3qYv^ + 



(l-4y2)g2 45 \ 33gy 5 



~ M 



{l+Y^)q^ 



1 + -v"^ -SqYv^ + 



{ -v^ - qYv^ + 



+ 2 qYv- 



(_7yi + iy 27\ 4 15qY 5 
^ K — 1 — V 



(i + y^y 9\ 4 



4^ =1 + ;^^^ + 



3 1 



(1 - y2)g2 27\ 4 ^ 5 

— H -3 qYv^ 



2 4 



1 + -i;2 + (2 - 3Y)qv'^ 



+eU^- -y + 2qYv'' + ( (1 - 2Y^)q^ - ) + 



99 



16 



43 gy 



^ (l + 7y)(l-y)92 ^ 27\ .4 ^ 3(2-5y)g ^., 



1 



+ + (4 - 7y)g?;^ 



.(3-15) 
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3.2. Reformulation of the source term for the Teukolsky equation 

Let us now discuss the method for rewriting the source term of the Teukolsky 
equation. The source term ()2-29p contains inverse powers of sin 9 in Cfhm, Cfhn-, and 
^-imifi ^ (sin0)~l™L These factors prevent us from performing the inverse Fourier 
transformation analytically. Therefore we want to remove all inverse powers of s'mO. 
Fortunately, the spheroidal harmonics Sa, C^Sa and C^C^Sa are proportional to 
(sine)!!""!-^!, (sin6l)ll™l-^l and (sin6l)l"'l, respectively. Therefore some of the inverse 
powers of sin^ cancel out, but we immediately find that some of them remain after 
this cancellation. In the following we show how we can eliminate all these annoying 
factors. 

For convenience, we introduce the new angular functions 



iS 



Sa 



(sin 0)1 



i5 = 



4Sa 



(sin 0)1^ 



^S = 



(sin6l)IH' 



and 



(3-16) 



fjSXri 



(sin ( 



\m\ ^—irmp 



X"", (m > 0) 

1, (m = 0) 
Xl"^!. (m<0) 



(3-17) 



If they are truncated at a finite post-Newtonian order, these new functions ^S, and 
hence a^m, have no inverse powers of sin 6, being expressed as polynomials of cos 9. 
There are a few special cases in which contains an additional overall factor of 
sind. We have 2S = sin^0(l it cos 9) x(a polynomial of cos 9) for m = ±1, while 
iS = sin^ X (polynomial of cos 9), and 2'S' = sin^ 9x (polynomial of cos 9) for m = 0. 
In addition, Xm can be expressed as a Fourier series of exp{if2g\), because it contains 
a positive power of X or X. Using a^m, the amplitude of the partial wave is rewritten 
as 



Za 



2iujB''^'' 



dXe 



-1 P 



2V27rp 



D2 

^ (o^m - 2iapi^rn) Ra 



2DrDg - (p + p)) + ia2SmiP " P)J~) Ra 



A 



+Dl2Sm{J--2pJ-)RA 



(3-18) 



r=r(A) 



where we have defined 



J- 
De 



dr - iK/A, 

£{r'^ + a^) -aC + 



aE 



C 



dr 
dX'- 



sm 



1 d cos 9 



sm 



dX 



(3-19) 



In this expression, it is clear that inverse powers of sin 9 are contained only through 
Dq, which contains the factor (sin^)^^, and therefore the index of the inverse powers 
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of sin 9 is at most four. We emphasize that this index is less for |m| = or 1, since 
2^171 and/or contain additional positive powers of sm9. Thus, we have to treat 
three cases, m = 0, m = ±1, and \m\ > 2, separately. Here we discuss only the 
case with \m\ > 2, and the other simpler cases with \m\ = or ±1 are deferred to 
Appendix A. 

Equation (13-18P has inverse powers of sin^ in DgiSm, Dg2^m, and Dl2^m- 
First we give a prescription for the terms containing the factor (sin^)~^, i.e. Dq iSm 
and Dg2^m- After some calculation, we obtain (for o" > 1) 



d „ dcosB „ 

dX dX 

d cos 6 
+ auj — h im 



dX 



where we have used 



C 



sm 



i ya£ 

C 

sin^ 

C 

sm 



C 



1 d cos ( 



sm a J sm 



dX 



+ 



d^W 
dX 



(3-20) 



(3-21) 



Using this relation, we can derive 



-^»n(o-— 1" 



m 1 a^m J i 



FraiA,B) = - 

m 



dcos 6 
dX 



dcos 9 



auj- 



dX 



im 



C 



sm 



(3-22) 



+ ima£ | B 



Since there are no inverse powers of sin ^ in FYni^a—i'^m^ cr^m)? the inverse Fourier 
transformation of this term can be easily performed. Obviously, the above formula 
is not valid for the case with m = because contains m~^. In (£/ sin^^)^ it 
seems that (sin0)~^ is remaining, but it is just a constant. 



C 



sm 



sm 



dcos 9 



1 



^e{cos9) 1 - cos2 9 



(3-23) 



which can be analytically integrated in the post-Newtonian expansion. 

Next, we consider the remaining term, Dg 2^m- Applying Eq. p-22p twice, we 
obtain 

DQ2^m = DeFmil^m, 2—m) 



FmiDg l^^rri, Dg 2"m) 7r-2^m 

m aX ^ 



m dX 



Using the geodesic equation for cos 9, one can show 

dDe 



dX 



cos 9{-D^ + 2ia£Dg + a^). 



(3-24) 



(3-25) 
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Substituting this relation into Eq. (I3-24P and applying Eq. (I3-22P again, we obtain 
cos 0^ 



m 

2ia£ cos 
m 



a COS U „ 



1 



m 



Sm)+^{1- Sin^ e)DI 2^m. 



Here, D| 2"m appears also in the last term on the right-hand side. After solving 
this equation for ^ we finally arrive at 



m 



2"m 



— 1 



1 + 



cos 9 



m 



Pm{Pm{o^mi 1 "m) ? -fm (l"m ) 2"m)) 



2ia£ COS 9 ^ , „ „ , cos ( 



m 



"-^m(l"m; 2"m) 



-2"m 



— 1 



(De sin^ 9)Fmil.^m, 2^m), 



(3-26) 



which is free from inverse powers of sinS. Note that this formula is not valid for 
m = and m = ±1 because of the factors of and (m^ — 

Combining the above results, the partial wave amplitude for \m\ > 2 is expressed 

as 



Za 



2iujB''^'' 



-1 -p 



2V27rP 



2Dr 



(o—m — 2iap\Em) Ra 



+ 



m? — 1 



{Fm{o^m, lSm){J- - [p + p)) + iaFmil^m, 2^m){p " P)J-} Ra 
{m + COs9)^mFmiFmioSm, lSm),FmilSm, 2^m)) 

- 2ia£ cos9 FmiiSm, 2^m) - cos9 2Sm^ 
-{De sin2 9)Fm{iS^, 2^™)! {J- - 2pJ-)Ra 



(3-27) 



r=r(A) 



This formula has no (sin0) ^ factor in the integrand. The formulas for |m| < 1 are 
presented in Appendix lAl 

§4. Results 

4.1. The evolution of orbital parameters 

Substituting Eq. (f3^ into Eqs. (fM2]) . ([2^ and (l2^ . the time-averaged 
rates of change for the three constants of motion up through 0{e^,v^) are given by 



. 73 2 

1 H 

24 



/ 1247 9181 ^2 I ,2 
V 336 672 
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dt 



'dC' 
. dt , 



+ 



73y 823y 



+ 



12 24 
/ 44711 172157 2 \ 4 



1375 



V 9072 ^ 2592 

329 527y2 r 4379 6533^2 ^ ^ 

"96 96~ ^ \T92 

/8191 44531 2^ _.5 



2 4 
q V 



+ 



V 672 ' 336 



(^)M.^(l-er^[( 



192 

/3749y 1759y 

TTV"' + 1 

V 336 56 

7Y 2^ /^1247y , 425y 2 \ ..2 



61 6iy^ 



+ 



63 



Y + 

9iy2 



V 336 



+ 



336 



24 8 [8 4 
/4471iy 302893y 2\ 4 



+ 



V 



V 9072 6048 

/57y 45y3 r20iy 37y3 



+ 



2 4 
7 V 



V 16 

/819iy 4836iy 2 
1 e I TTf" 

V 672 1344 

/2633 430iy2 f 66139 18419y2 , 2 



V 224 



224 



+ 



\ 1344 



448 



qv 



/743 23 2\ 2 /85y 21iy 
V336 ~42^ -|^ + 



+ 4 + 



+ 



97 



/ 129193 84035 2 \ 4 
V 18144 ^ 1728 ^ ^ 



/329 53y2 r 929 163y2 



V 96 8 1^ 96 
/2553y 553y n\ . 



2 \ 2 4 

e 1 g f 



V 224 



192 



4159 21229 
+ 



672 



1344 



(4-1) 



Here, a factor of (1 — e^)^/^ is factored out in order to make it easier to compare 
these results with the well-known formulas derived by Peters and Mathews. 

We can derive the evolution of the orbital parameters, = {u,e,y}, from the 
rates of change of the integrals of motion, /* = {f ,>C,C}, 



di' 
It 



h 



dP 
IF 



G] 



dP 



(4-2) 



Using the above relation with Eqs. (|3-4p and (|3-6|) . the evolutions of the parameters 
V, e and Y are obtained as 



' dv 
.'dt 



/743 55^2 
V336 21 
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/ 34103 _ 526955 2] 4 



V 18144 12096 
/329 815y2 r 929 

+ 1^ 



477y2 



de\ 
di/ 



V 96 96 
/1451F 1043y 2 

V 56 ^ " 



304 / M \ 8 



96 

;i_e2)3/2 



32 



z J / 

4159 48809 2^ 



/879y 515F 2\ 3 /^985 5969 2 



121 2\ /6849 4509 2^ 
~V2128"^^ " 



' 2128^ 



"76"' J"" ^\T52 
2064415 2^ „,4 



V 76 

/ 286397 2064415 2\ 4 

V 38304 ^ 51072 ' / ^ 

'3179 5869F2 J 8925 10747^^^ 2\ 2 4 
"608" 608 ^ \ 1216 1216 J ^ j ^ ^ 
'l903y 22373F 2\ 5 / 87947 4072433 2^ 



dY 
Itt 



8512 



V 304 

_244 /j. X 2 /2(1 . 

15 VM2; ^ ' ^ ' [\ 

fl3Y 277Y 2\ 2 4 / 10461 
~ V 244 244 ' ~ y 1708 



V 4256 68096 



1 + 



189 2\ 3 
e I gt; 



61 



^10461 83723 ^ 
3416 ' 



(4-3) 



Substituting Y = 1, it is found that these results are consistent with the previous 
results ^38)-40) gxcept for the errors reported in our previous paper. "^^^ 

4.2. Phase evolution of gravitational waves 

The time dependence of the orbital parameters gives us information about the 
phase evolution of gravitational waves. There are three fundamental frequencies, 
Ur = f2r {dt/dX)^^ , ujg = {dt / dX)^^ , and Lo^p = {d(p/dX) {dt/dX)^^. uJr is the 
angular velocity of radial oscillations, and it is not important for small eccentric 
orbits. By contrast, uiq and ui^ are equally fundamental for greatly inclined orbits. 
Stated briefly, gravitational waves cannot be expressed by a single phase, even if we 
pick up the dominant quadrupole contribution. 

When the orbit is in the equatorial plane, u,^ corresponds to the angular ve- 
locity used in the standard post-Newtonian calculation.^^) For general orbits, the 
post-Newtonian angular velocity is given by a; at = (1 — Y)uj0 + Yuj^. This can be 
understood in the following way. First we consider the frame rotating at the angular 
velocity of the orbital precession due to the spin-orbit coupling, uj0 — lj^. In this 
frame, the orbital plane of the geodesic motion is fixed. The orbital angular velocity 
observed in this frame is given by wg. Then, the effect of the rotating frame can 
be understood as the sum of two rotations, that around the axis perpendicular to 
the orbital plane and that around an axis in the orbital plane. The first rotation, 
whose angular velocity is Y{ojg — Utp), is understood as the relative motion between 
the rotating frame and the non-rotating one. Therefore, the angular velocity of the 
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orbit projected onto the momentary orbital plane is given by lon = {1 — Y)uj0 + Yujip. 
This is in some sense the angular velocity observed from the direction perpendicular 
to the orbital plane. The latter rotation causes the change of the orbital plane. As 
the orbital plane itself changes, this post-Newtonian angular velocity ujn does not 
appear in the waveform with a fixed direction of the observer. When we are con- 
cerned with a short period of time in which the precession of the orbital plane can 
be ignored, of course, there is no practical way to distinguish these three angular 
velocities, ujn, '^ip and uiq. 

In this subsection, we consider the evolution of the frequency F^p = lo^p/tt and the 
corresponding phase ^^{F^^) = 2ttJ Fy,dt, assuming that the evolution of all orbital 
parameters are governed by the time averaged expressions given in Eqs. ()4-3p . Here, 
we briefly describe the method for obtaining analytic expressions for them. The 
details, including the derivation of the Fourier transformed waveform under the 
stationary phase approximation and the results for Fg = ujq/tt, are reported in 
Appendix B. 

Using Eq. (j3-15p . we can express F^p in terms of . Solving this relation inversely 
with respect tov, we can express t; as a function F^, e and Y in the post-Newtonian 
expansion, i.e. expansion in powers of F^. Then, taking the derivative of F^{l^), we 
rewrite it using Eqs. (|4-3p . Substituting v{F^,e,Y) into this expression, we obtain 



dF^ 
dt 



10 73y 

\ h <! 

3 12 ^ 



1 H 

24 



743 3683 



336 ^ 112 



{ttMF^ 



{ 1193 
"36" 



+ I 4 + 
+ 



2335 



7r(7rMF^) + 



24 J 
/ 34103 



89353 



V 18144 1728 



(vrMF^)3 



233 527,^2 

VTY ^ Y^ + 

96 96 



2875 281,, 4165,., o 

\ Y H Y^ \ 

96 12 96 



qH^MF^ 



/4159 20135 



+ 



672 
743 
~72 



+ 



96 
13907y 

336 



7r{7rMF^y- 



+ 



16687 84365y 1 



84 



336 



e2 qiirMF^ 



(4-4) 



where we have introduced the chirp mass, M = M^/^/i'^/^. The evolution equation 
for e and Y can also be rewritten by using F^p instead of v. We have to integrate 
these three differential equations simultaneously. This may appear to be difficult, 
but we can perform this integration iteratively, assuming a small eccentricity and 
the post-Newtonian expansion. It is important to note that 5Y = Y — Yj is also 
small in the post-Newtonian expansion, where Yj is the value of Y at F^ = 0. For 
this reason, in the lowest order approximation, we can fix the value of Y to Yj. The 
effect of the variation 5Y can be taken into account perturbatively. Eliminating dt 
from these differential equations, and Y are integrated as functions of F^p as 



= ejiirMF^ 



19 
9 



l + ^(^^MF^] 
1008^ ^' 
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377 /38 287,, , , , 



+ 



r 35705555 9 /731 19,, 1675, . ^.^x* 

< h g Yj H Y? } (7rMF^ )3 

\ 1524096 ^ \576 6 576 V i 

269971 /41795 5435639, ,,^,5 . 



= _(1 _ y/)<^(vrMF^) ( ^ - ige1(vrMF^)-f +•••). (4-5) 



Here, e/ denotes the limiting value of e x (ttMF^)^^/^^ in the limit 0. We have 

denoted the leading-order correction to due to the evolution of Y separately by 
6e'^^Y)j although it is relatively order 3PN. 

Substituting the relations in Eq. (j4-5p into Eq. (j4-4|) and integrating F^{dFy,/dt)^^ 
over F^, the phase 'P,p{F^) is calculated as 



+ 



^i^MF^yl 

/25 365,, 
15293365 o 



l + ^(^MF 



1016064 
r 38645 
\ 672 



1008 
- lOTrj {ttMF^] 
1165 2635, 



WYi 

96 

/3715 688405 
V 168 2016 



y/ \ (7rMF^)3 ln(7rMF^ 



/ 22265 



864 



ejiirMF^)- 

r 65561 

1 2880 
111064865 

10948608 
r 3873451 



785 
'272 



21072205 , 
2045665 



(vrMF, 



225792 



{ttMF^ 



/2057 2953,, 
+ ^^1/ 



+ 



90 540 

698695 16895 



(ttMF^) 

650665, 



\ 86184 



-vr 



101376 
/ 1247185 



-Yi + 



24624 



+ 



1056 101376 
1899015067 



8273664 



y/ ^ {ttMF^ 



(4-6) 



where is a constant of integration, and 5^^^ and (5<?^^^^ express the corrections 
of 0{ej) and the terms associated with the time variation of the inclination angle, 

respectively. We have included the cross terms between these two effects in (5^^^^ . 

Here, several remarks are in order. The terms of 0(e|) seem to be large in 
the sense of the post-Newtonian order. Those terms have relatively large inverse 
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powers of F^. However, we should recall that ~ ej{7rMF0)~^^^^ . Therefore, 
under the current assumption of a small eccentricity (e <C 1), the terms associated 
with the factor e'j{TrMFg)~^^^^ are much smaller than the other terms. In the above 
calculation, we have basically kept the terms up through 0{v^,e'^). However, for 

6^lp'^\ we have also kept the terms of higher post-Newtonian order, 0{v^), since the 
leading-order correction due to 5Y starts with this order. 

As mentioned above, more explicit formulas including the Fourier transform of 
the waveform under the stationary phase approximation are given in Appendix B. 
Here, we would like to stress that the analytic method for computing the phase 
evolution presented here can be extended systematically to higher order in v and e. 

§5. Summary 

In this paper, we have improved an analytic method for calculating the adiabatic 
evolution of the orbital parameters of a point particle orbiting a Kerr black hole. We 
have removed the previous limitation to an orbit with a small inclination angle. To 
do this, we rewrote the source term of the Teukolsky equation, removing inverse 
powers of sinOz by using integration by parts, where 9z is the polar angle of the 
particle position. As a result, our new expression for the source term is much more 
complicated than the original one, but it is possible to carry out an analytic Fourier 
transform of this source term. Moreover, even in a numerical approach along the 
line of Refs. 42) and 43), our expression might be useful when we consider the polar 
orbit, since the singular behavior of the source term near the poles in the original 
formulation is completely removed. 

We have also shown that it is possible to integrate the time evolution of the 
orbital parameters under the assumption that it is governed by the time averaged 
values derived in the adiabatic approximation. Furthermore, analytic calculation of 
the Fourier transform of the waveform is also demonstrated in Appendix B. The 
results presented here are restricted to 0{v^,e'^), where v and e are, respectively, the 
velocity and the eccentricity of the particle. Our method, however, can be extended 
systematically to higher order. Such analytic expressions for the waveform should 
be useful for the fast generation of templates of gravitational waveforms. 

Our analytic expressions for the averaged evolutions of orbital parameters ()4-3p 
are written as polynomial functions with respect to Y defined by (|3-6p . This result 
could be anticipated from the expressions for orbits given in ^3.11 The only non- 
polynomial expression is Vl — in cos 0(A). However, changing the signature of 
cos 0(A) alone corresponds to simply flipping the direction of the z-axis. This change 
should not affect the rates of change of the energy, the angular momentum and 
the Carter constant, which means that only even powers of cos 0(A) contribute to 
these rates of change. Here we do not pursue a rigorous proof of this statement. If 
we can prove that the averaged evolutions of orbital parameters can be written as 
polynomials in Y, we will not have to perform computations for a large inclination 
angle, at least for analytic calculations, because the expressions obtained in the 
expansion of a small inclination angle will be sufficient to determine all the coefficients 
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that we want to know. Then, extrapolation to a large inclination angle is exact. To 
do this, we also need to know the maximum power of Y at each post-Newtonian 
order. This issue will be studied in a future publication. 
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Appendix A 

Reformulation of the Teukolsky equation for |m| < 1 



m = 0. In this case, (finite polynomial of cos 9) in the PN 

approximation. Thus, we can write Eq. (j3-18p as 



Za 



2DriD0 sm^9) 



CO 

2 , 



-1 p 



2V27rp 



[oSm - 2iapi~rn) Ra 



+ {Do sin 



A \sm 

2 n\2 , ^2 



sm 



sm 



{Jl - 2pJ^)Ra 



(Al) 



r=r(A) 



This form has no factors of 1/ sm.9. 
• m = ±1. In this case, the DgiSm and D'^2^m terms have l/s\v?9 factors. 
We can treat the Dq iSm term here in a manner similar to that for \m\ > 2. 
Therefore, we have to consider the Dg 2^m term. Only in the case m = ±1, we 
can find the equation 



dX^^' 



De + 



1 



dcos ( 



ia8 + i 



C 



1 ± cos 9 dX \ sm n /^^ 

Since 2»S' is proportional to sin^ ^(1 it cos 9), we obtain the relation 



A±i. (A-2) 



Di 



{De sin^ 



2S 



sm''9{l±cos9) 



+(1 ±cos6') ±— + ia£ 
dX 



d 



dcos t 
dX 
C 
sm 



A±i. (A-3) 
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Therefore, the amphtude of the partial wave for m = ±1 is 



d\ e 



-1 p 



(o—m — 'iiapi^m) Ra 



2V27rp, 

Fmio^m, l^m){J- - {p + P)) 



+ia{Dg s\v? 



2—m 



sm 



{-p-p)J-\RA 



+ 



{De sin^ 9) 2S 
sin^ 9{1 + m cos 9) 



dcos 9 
dX 



Xr, 



+(1 + m cos 9) m— + ia£ 
dX 



d 



sm 



^[JI-2pJJ)Ra 



r=r(A) 



(A-4) 



Appendix B 

Formulas for phase evolution 



In this appendix we present all the detailed formulas for the orbital and phase 
evolutions to supplement §4.2[ The expression for dFo /dt corresponding to (j4-4p is 



dt 



96 , ..^sii 
:{7rMFg)ir 



f73Y , 857Y 2 



, 157 n\ /743 3683 



H e- 

V 336 112 



{■KMFg) 



12 



+ 



24 



e' q{TTMFe) 



/ 34103 89353 



+ 



48 

233 527, o 

\ + 

96 96 

/4159 20135 



V 18144 1728 



2875 4165, 



96 



+ 



96 



Y^\e^]q^{TiMFeY 



672 96 ' ^ 



/13907y 84365y 



+ ) q(-KMF0)t 

336 336 ' ^ ^ 



(B-l) 



As mentioned in ^4.21 the standard post-Newtonian calculation uses a frequency that 
is defined differetly, namely Fjv = (1 — Y)Fq + YF^. For comparison, we also give 
an expression for dFj^/dt, 



dFN 
dt 



57r7W2 
f 34103 

+ 



l-^(7rMF^)i + (47r 



ii3y 
12 



(vrMFiv) 



\ 18144 



233 



527, 



+ -— +2Y + —Y' \ {ttMFn) 



96 



96 



(B-2) 



This agrees with the standard post-Newtonian result, up to 1.5PN.^^) 
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The expression for the phase evolution of the ^-osciUation corresponding to (j4-6p 



IS 



(SY) 



= ejinMFe) 



^ + lui^^^^^)' + {^^^^ - ^ ^^^^^ 
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+ 



1165 2635 



1016064 " V 96 
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225792 
Yi i (vrMFg) 



(TiMFe 



+ 



10948608 ■ V101376 101376 
r 3873451 1899015067 



\ 86184 



8273664 



Yi } (TTMFe)- 



(B-3) 



where (pg is a constant of integration. 

In order to compute the waveform, we also need to know t{F) = J {dF/dt)^^dF. 
This computation can be carried out in a manenr similar to that for the calculation 
of the phase ^{F). The results are 
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where and are constants of integration. 

Writing the Fourier components of the waveforms as h{f) = Ae'^'^^f\ V'(/) is 
calculated as ip{f) = — #(/) + 27r/i(/), using the stationary phase approximation.^^) 
The final results are 
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Here, we define 6t^ = tg — t^, and we have absorbed a constant phase into (p^ and 
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